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Rotor blades in a radial turbine with nozzle guide vanes typically experience harmonic aerodynamic excitations 


due to the rotor stator interaction. Dynamic stresses induced by the harmonic excitations can result in high cycle 


fatigue (HCF) of the blades. A reliable prediction method for forced response issue is essential to avoid the HCF 


problem. In this work, the forced response mechanisms were investigated based on a fluid structure interaction 


(FSI) method. Aerodynamic excitations were obtained by three-dimensional unsteady computational fluid dy- 


namics (CFD) simulation with phase shifted periodic boundary conditions. The first two harmonic pressures were 


determined as the primary components of the excitation and applied to finite element (FE) model to conduct the 


computational structural dynamics (CSD) simulation. The computed results from the harmonic forced response 


analysis show good agreement with the predictions of Singh's advanced frequency evaluation (SAFE) diagram. 


Moreover, the mode superposition method used in FE simulation offers an efficient way to provide quantitative 


assessments of mode response levels and resonant strength. 
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Introduction 


Nowadays, turbocharger plays an important role to 
improve fuel economy in internal combustion engines. In 
order to improve operational range and performance effi- 
ciency of the turbocharger, nozzle guide vanes are often 
employed to control the flow of exhaust gas to the radial 
turbine. However, using of the stator vanes is inevitable 
accomplished with the harmonic aerodynamic excitations 
caused by the interaction between the stator vanes and 
rotor blades. The harmonic excitations can excite natural 
frequencies within the turbine rotor, resulting in reson- 
ance of the rotor blades. Dynamic stresses induced by the 
harmonic excitations can lead to HCF and eventually 
result in failure of part [1, 2]. To ensure the operation 
safety of turbocharger, the forced response issue has be- 
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come the main concern for the reliabilities of radial tur- 
bines. 

Avoiding resonance conditions is one of the most crit- 
ical issues regarding bladed disk design. The most simple 
and common approach to prevent resonance excitations 
is the Campbell diagram criterion. However, this ap- 
proach is too conservative for bladed disk system, which 
has modes with very close natural frequencies [3]. A 
more effective criterion for bladed wheel design is re- 
presented by SAFE diagram. The importance of both 
frequency and shape matching conditions is explained in 
the book written by Bloch and Singh [4]. The SAFE dia- 
gram method was applied in more recent papers and its 
effectiveness has been demonstrated [5, 6]. 

The SAFE diagram method is a qualitative analysis 
approach; thus, it cannot provide information about re- 
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Nomenclature 
An Fourier coefficient 
[C] damping matrix 
dn ND of the n-th harmonic pressure 
{F} harmonic load vector 
F, real load component 
F> imaginary load component 
[K] stiffness matrix 
[M] mass matrix 
m,n harmonic index 
Ng number of rotor blade 
Ny number of stator vane 
Pio time averaged value of pressure (Pa) 
Pin real pressure component of the n-th harmonic 
P>, imaginary pressure component of the n-th 
; harmonic 
Pr pitch length of rotor (m) 
Ps pitch length of stator (m) 
t time (s) 
T period of one pressure fluctuation cycle (s) 
At time interval (s) 


sonant strength. To evaluate the dynamic stresses or 
strains in bladed disk system, utilizing experimental test 
method is a credible option. However, it has distinct 
drawbacks in terms of test preparation time and cost to 
conducting an experiment (e.g., by applying strain gaug- 
es or using tip-timing methods). The FSI method, an in- 
terdisciplinary analytical approach which coupled CFD 
and CSD analysis, has been developed to speed design 
cycles. Besides the time and cost advantages, transient 
CFD simulations provide opportunity to investigate the 
character of aerodynamic excitations causing resonance. 
Moreover, the stresses and vibration amplitudes of struc- 
ture can be derived from CSD analysis. The FSI method 
has been presented and validated in several research 
works [7-10]. 

To achieve improvements of resonance prediction 
quality and avoid HCF problems early in the design pro- 
cess, further knowledge of the underlying forced response 
mechanisms that govern HCF failures is still desired. In 
this paper, a weakly coupled FSI analysis was conducted 
to investigate the forced response behavior of radial tur- 
bine. Aerodynamic excitations were obtained by a three- 
dimensional unsteady CFD simulation. Harmonic pressure 
data including the frequency, amplitude and phase angle 
information was transferred to the FE model to perform 
structural dynamics analysis. The results were analyzed 
and compared with the predictions of SAFE diagram. 
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U; flow condition at the lower boundary 
Uy flow condition at the upper boundary 
Vr velocity of rotor (m/s) 

{x} nodal displacement vector 

Greek symbols 

y phase angle of load (degree) 

@ circular frequency of rotor (rad/s) 

Q excitation frequency (Hz) 

Q' fundamental frequency (Hz) 
Abbreviations 

CFD Computational Fluid Dynamics 

CSD Computational Structural Dynamics 
EO Engine Order 

FE Finite Element 

FSI Fluid Structure Interaction 

HCF High Cycle Fatigue 

ND Nodal Diameter 

SAFE Singh’s Advanced Frequency Evaluation diagram 
VPF Vane Passing Frequency 


Radial Turbine Description 


The radial turbine under investigation was designed to 
achieve an expansion ratio of 1.75 at rotational speed of 
63200 r/min. The main geometrical parameters of the 
radial turbine are summarized in Table 1. The material 
used for the turbine rotor is Fe-Ni-Cr alloy, the density is 
7860 kg/m’, the Poisson’s ratio is 0.3, the Young’s mod- 
ulus is 1.1<10° MPa, and the yield strength is 726 MPa. 


Table 1 Geometrical parameters of radial turbine 


Impeller blade number 13 
Stator vane number 15 
Impeller inlet diameter (mm) 96.5 
Impeller tip width (mm) 14 
Impeller exit hub diameter (mm) 21 
Impeller exit shroud diameter (mm) 77 
Inner diameter of nozzle (mm) 100 
Outer diameter of nozzle (mm) 136 


CFD Model 


In this study, the unsteady flow field information was 
obtained by transient CFD simulation, which was im- 
plemented in a time-implicit Navier-Stokes equation 
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solver in combination with shear stress transport turbu- 
lence model. The radial turbine has prime number of 
blades. To get accurate results, it need to model all the 
turbine passages to reduce the pitch scale ratio to unity. 
This would be a very expensive calculation. For this rea- 
son, the rotational phase shifted periodic boundaries were 
used to enable only a small section of the full geometry 
to be modeled. The phase-shifted periodic condition is 
based on the principle that the flow condition at the upper 
boundary of the stator Uy at a certain time tot+At is the 
same as at the lower boundary of the stator U; at time 
instant tọ [11]. The periodicity of the boundary condition 
in phase shift is illustrated in Fig. 1 and can be expressed 
in the following equation [12]: 
U; (ty) =Uy (ty + At) (1) 
The time interval At that the lower boundary lags the 
upper boundary is defined by (Pr-Ps)/Vp. Here Pp and Ps 
are the pitch lengths of the rotor and stator, respectively, 
and Vp is the velocity of the rotor. 


Stator Rotor Stator Rotor 


Uv V; Pr Uv 
“om |p EBs 


t=to t=to+At 


Fig. 1 Phase shifted periodic boundary condition 


Based on the phase shifted periodic boundary condi- 
tion, the Fourier transformation method was applied for 
efficient calculation of the unsteady aerodynamic flow. 
The Fourier transformation method makes use of tem- 
poral Fourier series decomposition (Equation 2) at the 
fundamental frequency @to avoid storing the flow signal 
for a full time period on all pitch-wise boundaries, simi- 
lar to the method proposed by He L [13]. 


M , 
f (t) 2 X Age) (2) 
m=—M 


Where m is the harmonic index, æ is circular frequen- 
cy of the turbine rotor, A,, is the Fourier coefficients. 

Only the Fourier coefficients A,, need to be stored to 
reconstruct the flow solution at an arbitrary time and 
pitch-wise location. The Fourier transformation method 
uses a double-passages model, where the Fourier coeffi- 
cients are sampled and collected at the middle interface 
between the two adjacent passages. The fluid mesh used 
for the double-passages approach is shown in Fig. 2. 


CSD Model 


To determine vibration characteristics of elastic struc- 


J. Therm. Sci., Vol.25, No.2, 2016 


tures, a modal analysis can be performed first. The linear 
equation of motion for free vibration is written as: 
[M ]{x} +[C] {x} +[K] {x} = {0} (3) 

Where [M] is structural mass matrix, [C] is structural 
damping matrix, [K] is structural stiffness matrix, {x} is 
nodal displacement vector. 

The response of the structure under harmonic loads is 
obtained by solving the general equation: 

[M ]{x} +[C] ix} +[K] {x} = {F} (4) 

Where {F} is applied harmonic load vector. 

By using Euler's identities, the harmonic load applied 
on the structure with different phases can be represented 
compactly with complex notation: 

{F} = {Fpe" je = {F + jF }e™ (5) 

Where Fm is magnitude of the load, y is phase angle 
of the load, Q is excitation frequency, F, is real load 
component and F; is imaginary load component. 

Based on the assumption that the structure vibrates 
with harmonic motion, the governing equation (3) can be 
converted to an eigenvalue equation. The Block Lanczos 
method was applied for searching eigenvalues (natural 
frequencies) and eigenvectors (mode shapes) of the equa- 
tion. Then the mode superposition method was employed 
to solve the general equation (4) in modal coordinates. 
The dynamic response of a structure to harmonic excita- 
tion can be characterized by using the eigenvalues and 
eigenvectors obtained from the modal analysis. In mode 
superposition method the nodal displacement is calcu- 
lated from a superposition of eigenvector [14]. 


Sampling Interface 


Fig. 2 Double-passages fluid mesh 


For structures with cyclic symmetric nature, structural 
dynamics analysis can be performed for the entire struc- 
ture by modelling only one sector of the structure. The 
basic sector is duplicated in the dynamic analysis to sa- 
tisfy the cyclic constraint relationships and to obtain 
nodal displacements. This technique was introduced by 
Dickens [15]. Fig. 3 shows the FE mesh of the basic 
radial turbine rotor sector used in this study. 
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Fig. 3 FE mesh of the basic rotor sector. 


Unsteady Aerodynamic Loads 


The mechanism responsible for the unsteadiness in the 
flow field was studied by transient CFD simulations. As 
regards this radial turbine, a rotor blade interacts 15 times 
with the stator vanes after a complete revolution of the 
impeller; thus, a 15th engine order (EO) excitation can be 
expected. EO represents the number of pulses per-revolu- 
tion. A representative position at mid-span was chosen to 
investigate the unsteady behavior of flow within the tur- 
bine passages. Fig. 4 shows the entropy contours of flow 
field at four time points covering the time interval of ro- 
tor blade passing. Here T is the time a blade needed to 
cross one pitch length of stator. The stator wake imping- 
ing on the blade surfaces is a significant source of the 
periodic pressure fluctuation in the rotor blade passage. It 
can be observed that the stator wake is chopped into short 
segment by the moving rotor blade. The chopped-off 
wake segment migrates toward the suction side of the 
blade and decays downstream. After being cut by the 
stator vane, the chopped wake segment drifts across the 
blade passage due to the negative jet effect. The chopped 
wake segment acts as a negative jet because the tangen- 
tial velocity component of the wake flow is higher than 
the mean flow in reference frame of rotor [16]. The 
chopped wake segment migrates with direction from the 
suction side to the pressure side and subsequently mixes 
with the vortex flow near leading edge of the blade. The 
unsteady velocity field induced by the stator wakes 
strongly affects the pressure field of blade passages. The 
contours of fluctuating pressures at the same four time 
points are shown in Fig. 5. Another mechanism which 
significantly influence the pressure field of blade passag- 
es is the potential field interaction. The potential fields of 
isolated stator row and blade row are steady in their re- 
spective reference frame. Since the blade row has a rela- 
tive motion with reference to the stator row, the spatially 
non-uniform potential field of stator row can also cause 
the pressure fluctuation in the blade passages. In order to 
study the variation of the pressures from one time to the 
next, three monitoring points (P1, P2 and P3, as shown in 
Fig.5) were selected along the blade chord direction. 
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Fig. 4 Entropy contours at mid-span 
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Fig. 5 Pressure fluctuation and the positon of monitoring 
points 


Fig. 6 shows the curves of the pressures changed with 
time at the monitoring points. It can be found that the 
variation curves exhibit periodic changing characteristics. 
These variation pressures have the same period, but they 
are out of phase with each other and have different am- 
plitude values. Monitoring point P1 located near the 
blade leading edge has the highest averaged pressure, 
followed by the monitoring point P2, the monitoring 
point P3 has the minimum averaged pressure. The aver- 
aged pressure decreases along the blade passages because 
the gas expands adiabatically and the pressure energy is 
transferred to kinetic energy. 

The pressure variations at the fluid mesh nodes of 
blade surfaces are periodic functions of time t; therefore, 
these functions can be expanded in Fourier series. By 
using complex notation to implicitly consider both am- 
plitudes and phase angles of the pressure loads, the Fou- 
rier series of pressure functions at blade surface nodes 
can be represented as: 


N Ta 
PO =R +D ((An+ ine") © 
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Fig.6 Variation of the pressures 


Where n is the harmonic index, Pı n and P2,, is the real 
and imaginary pressure components of the n-th harmonic, 
Q’ is the fundamental frequency, and P; is the time av- 
eraged value. 

Since the period of one pressure fluctuation cycle is 
the time a blade needed to cross one pitch length of stator, 
the fundamental frequency of Fouries series (6) is identi- 
fied equal to the vane passing frequency (VPF). VPF can 
be calculated by the circular frequency of turbine rotor 
times the number of stator vanes. 

To perform CSD analysis subsequently, the pressure 
data at fluid mesh nodes was transferred to solid finite 
elements by proximal interpolation algorithm. The dis- 
tributions of the first two harmonic pressures at blade 
surface elements are shown in Fig. 7 and Fig. 8, respec- 
tively. According to the Fourier decomposition, the first 
harmonic pressures (Pı ı and Pz; at 15th EO) and the 
second harmonic pressures (P2 and P22 at 30th EO) 
have much larger amplitudes than the high-order (n 2 3) 
harmonic pressures; thus, the first two harmonic pres- 
sures were determined as the primary components of the 
aerodynamic excitation. 


Modal Analysis 


The eigenvalues can be changed by the stress stiffen- 
ing because of the fact that the stress state would influ- 
ence the values of the stiffness matrix. Therefore a pre- 
stressed modal analysis was conducted on the turbine 
rotor considering the initial stress state, which included 
the centrifugal stress caused by rotor spinning and the 
stress caused by time averaged aerodynamic pressure. 
Fig. 9 shows the equivalent stress distribution of the rotor 
under these two static loads. The maximum static stress 
occurs at roots of the blades and shows compliance with 
strength requirements. 

Mode shapes of cyclic symmetry structures are cha- 
racterized by specifying the number of nodal diameters 
(ND). The SAFE diagram can be used to determine the 
critical frequency and corresponding ND. The basic idea 
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Fig. 7 Distribution of the first harmonic pressures 
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Fig. 8 Distribution of the second harmonic pressures 
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Fig.9 Static stress distribution 


behind SAFE diagram is that the resonance happens 
when the frequency and periodic shape of harmonic ex- 
citation match the natural frequency and mode shape of 
the structure. The shape matching conditions have been 
analytically derived in reference [3]. The matching num- 
ber of ND for the n-th harmonic excitation can be calcu- 
lated by the following equation: 


d, = ny - Np| (nNy +| N; 12 J) J/N] (7) 


Where d, is ND of the n-th harmonic excitation, Ng is 
the number of rotor blade, Ny is the number of stator 
vane, “| J" is the notation of the floor function which 


find the integer part of real number. 

The ND of the first two harmonic pressures was cal- 
culated by using equation 7. The high-order harmonic 
pressures were neglected in this study because the am- 
plitude is small. Regarding the radial turbine with Ng = 
13 and Ny = 15, the results dı = 2 and d, = 4 can be ob- 
tained. 

In order to show a clear view of the matching condi- 
tions, a SAFE diagram (as shown in Fig. 10) was con- 
structed, which comprises the information of natural fre- 
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quency, ND and the first two order harmonic pressures. 
The solid points in the SAFE diagram represent the nat- 
ural modes arranged by frequency and ND. Because the 
maximum number of ND is not larger than half of the 
blade number , this turbine rotor with 13 blades has up to 
6 NDs. The excitation line crossing frequency of each 
EO excitation was also plotted with dashed line. Since 
the CFD analysis was conducted for the design speed 
condition, the turbine rotor running at 63200 r/min has 
the first harmonic pressure occurring at 15800 Hz, and 
the corresponding EO is 15. The second harmonic pres- 
sure occurs at 31600 Hz, and the corresponding EO is 30. 

The SAFE diagram shows that there is no full reson- 
ance excited since no frequency is coincident with the 
first harmonic pressure or the second harmonic pressure. 
However, mode 6 and mode 7 of 2 ND modes have fre- 
quency value close to the first harmonic pressure and the 
shape matching is verified. Thus the excitations of these 
two modes still can be expected. Similarly, the mode 15 
and mode 16 of 4 ND modes can also be verified as the 
critical modes for they have the same shape and close 
frequency value with the second harmonic pressure. 

Fig. 11 shows the two critical 2 ND modes obtained 
from modal analysis. The appearances of the turbine ro- 
tor all contain areas of zero displacement where the nodal 
lines cross the entire disk. Although these two modes can 
be classified as 2 ND modes and have the same mode 
shape, the deformation forms of rotor blades are different. 
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Fig. 10 SAFE diagram 
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Fig. 11 Critical modes 
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The phenomenon that local regions of the impeller exhi- 
bit dominant mode shapes and frequencies nearly iden- 
tical to that of a blade fixed alone has described in the 
reference [17]. Thus a particular impeller mode is consi- 
dered to be the modal coupling of the impeller back wall 
and blades. 


Harmonic Forced Response Analysis 


To have an insight into the dynamic behavior of the 
rotor blade system, a harmonic response analysis was 
carried out by utilizing the mode superposition method, 
which provides an effective solution for multi degree-of- 
freedom system dynamics problems. In mode superposi- 
tion method, total response of a system can be obtained 
by summation of the mode shapes, each multiplied by a 
modal coefficients. The values of modal coefficients re- 
present the response amplitudes applied to each mode shape. 

The aerodynamic loads obtained from the unsteady 
CFD results were applied as boundary conditions of the 
finite element solutions. According to reference [7], aero- 
dynamic damping was ignored because of the small blade 
size. The modal coefficient histogram of the modes ex- 
cited by the first harmonic pressure was plotted in Fig. 12. 
The plot was limited only to display the modal coeffi- 
cients of 2 ND modes because only these modes have 
response to the first harmonic pressure. The results show 
that mode 6 and mode 7 are significantly excited, which 
is consistent with the prediction of SAFE diagram. Fur- 
thermore, the influence of each mode on the modal re- 
sponse can be evaluated. Although mode 6 and mode 7 
play more prominent roles in the overall response, the 
lower order modes are also responsive to the first har- 
monic pressure. However, for modes with high natural 
frequency, the responses are small. For example, mode 2 
and mode 3 are excited considerably but there are less 
responses for the mode order larger than 12. 

Fig. 13 shows the modal coefficient histogram of 4 
ND modes excited by the second harmonic pressure. It is 
found that the response behavior has the similar character 
compared with that of the first harmonic pressure. Only 
the 4 ND modes have response to the second harmonic 
pressure due to the shape matching conditions. Likewise, 
mode 15 and mode 16 are significantly excited, corres- 
ponding to the expectation of SAFE diagram. In general, 
the contribution of the modes which have low frequency 
or close frequency value with the excitation are more 
pronounced than the others. 

The spatial distribution of the harmonic pressures is 
another key factor affecting the response of the modes. 
The effect is reflected in that some modes (e.g., mode 4, 
mode 7 and mode 11 in Fig.13) have relatively small 
response to the excitation. As mentioned above, the blade 
deformation form of a particular mode is different from 
others; therefore, the distribution of the pressure loads on 
blade surfaces should be considered. In most cases, it is 
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hard to estimate the response level by physical intuition. 
The mode superposition method provides a useful ap- 
proach to estimate which modes will be excited effective- 
ly by the spatial distribution of the harmonic pressures. 
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Fig. 13 Modal coefficients of 4 ND Modes 


Although the SAFE diagram is able to indicate the 
critical modes that encounter forced response, it provides 
no information about the resonant strength of structure. 
The application of mode superposition method in finite 
element analysis not only can judge whether the reson- 
ance conditions are acceptable, but also obtain the results 
of resonant amplitudes and dynamic stresses in detail. 
The contour of dynamic stress induced by the first har- 
monic pressure is shown in Fig. 14. It can be seen that 
the dynamic stress distribution is in a pattern with 2 NDs 
presenting on the disk because the resonance is contri- 
buted by the excited 2 ND modes, and the maximum 
stress is located on the leading edges of the rotor blades. 
Fig. 15 shows the contour of dynamic stress induced by 
the second harmonic pressure. The dynamic stress distri- 
bution is in a pattern with 4 NDs presenting on the disk 
and the maximum stress is also located on the leading 
edges of the rotor blades. The dynamic stress induced by 
the first harmonic pressure is six times larger than that by 
the second harmonic pressure, which could be considered 
as the major alternating stress on the rotor. 
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Fig. 14 Dynamic stress distribution under the first harmonic 
pressure 
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Conclusions 


A FSI method adapted for the forced response analysis 
of bladed disk subjected to excitations has been applied 
to a radial turbine. The major goal of the approach is to 
quantify the forced response of the turbine rotor under 
aerodynamic excitations caused by the rotor stator inte- 
raction. From the analysis results, following conclusions 
can be derived. 

The periodic pressure fluctuations in the rotor blade 
passages are arising from the wake interaction and the 
potential field interaction. The periodic functions of sur- 
face pressure can be expanded in Fourier series with the 
fundamental frequency equal to VPF. According to the 
Fourier decomposition, the first two harmonic pressures 
were determined as the primary components of the exci- 
tation. 

A SAFE diagram was constructed to predict critical 
modes of the turbine rotor. Mode 6 and mode 7 of 2 ND 
modes are verified as the critical modes because the fre- 
quency values are close to the first harmonic pressure 
and the shape matching is verified. And mode 15 and 
mode 16 of 4 ND modes are likewise verified as the crit- 
ical modes. 

The computed results from the harmonic forced re- 
sponse analysis show good agreement with the prediction 
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of SAFE diagram. The critical modes are significantly 
excited and play a more prominent role in the overall 
response. Quantitative assessment of the mode response 
levels was obtained. On the whole, the contribution of 
the modes which have low frequency or close frequency 
value with the excitation are more pronounced than the 
others. 

It is found that the maximum dynamic stresses in- 
duced by the first two harmonic pressures both are lo- 
cated on the leading edges of the rotor blades. The dy- 
namic stress induced by the first harmonic pressure is 
considered as the major alternating stress. These stress 
data will be used to perform fatigue and durability analy- 
sis in the next stage of research work. 
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